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Justin R. Mason 
Quantum Corrections to Diffusion in Stars 



Quantum corrections can be important for diffusion and the melting temperature 
of dense plasmas in compact astrophysical objects, particulary white dwarfs and 
neutron stars. Typically ions in these systems are modeled classically, but Daligault 
et al. use a semiclassical inter-ion potential. We run molecular dynamic simulations 
using this semiclassical approach in order to calculate the diffusion coefficient and 
melting temperatures in a one component plasma. We find that in liquid simulations 
quantum corrections do not have a significant effect on diffusion, increasing it by only 
a factor of two from our classical simulation. However, in solid simulations, diffusion 
slowly increases for small quantum corrections, but once quantum effects become 
large enough, the system can liquify. We also find that quantum corrections can 
decrease the melting temperature of a one component plasma, potentially affecting 
the way in which white dwarfs cool. These results suggest that a helium white dwarf 
may remain liquid at typical white dwarf densities, while the structure of a neutron 
star's crust could be altered due to quantum effects. 
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Chapter 1 
Introduction 



By understanding quantum effects in diffusion in white dwarfs (WD) and neutron 
stars (NS), one can determine more accurately the structure, age, and evolution 
of the Milky Way Galaxy. In a WD, sedimentation of ions with a large mass-to- 
charge ratio releases gravitational energy and slows their cooling rates [l|. The crust 
of a NS undergoes chemical separation during crystallization at extreme densities 
This separation changes physical properties of the crust. Diffusion constants 
for a system can be determined from Molecular Dynamics (MD) simulations. As 
computational power has increased a lot of effort has been put into MD simulations to 
determine diffusion constants for these Coulomb plasma systems. These simulations 
are extremely powerful tools that can be used to understand the underlying physics 
for systems such as WDs or NSs by matching observational evidence. 

WD stars are the remnants of low- and intermediate-mass stars that can no 
longer sustain nuclear fusion and are believed to mark the end of stellar evolution 
for more than 97% of stars [jij]. WDs spend their lives slowly cooling by radiating 
thermal energy from their cores Qj. WDs are extremely compact and dense stellar 
objects. Their mass is on the order of that of the Sun, but they are approximately 
the size of the Earth. Reaching n ~10 28 -10 30 cm -3 densities, they are comprised of a 
pressure ionized Coulomb plasma, which is stopped from further collapse by electron 
degeneracy pressure. In such a degenerate gas the electrons are distributed fairly 
uniformly surrounding the nulcei. However, due to strong Coulomb interactions at 
high density the ions undergo a phase transition to form a plasma crystal. Any 
energy released due to gravitational collapse is used in raising the Fermi energy 
of the electrons (forcing degenerate electrons into higher energy levels) rather than 
increasing the star's luminosity Q. 

The WD luminosity function is the number of WD stars observed as a function of 
magnitude. Structure within the WD luminosity function (LF) provides information 
on the age and timescale of formation for components of the Milky Way Galaxy. In- 
terpretations of the WD LF structure is dependent upon our underlying knowledge 
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of cooling physics (3|. In the field of WD cosmochronology, these slowly changing 
stars are used to track the age and star formation history of the Milky Way. A good 
review of WD cosmochoronology is given by Fontaine et al. Precise determi- 
nations of the WD LF in globular cluster NGC 6397 has shown a peak around a 
magnitude of 26.5. A build up of WDs at this magnitude is attributed to the release 
of latent heat of crystallization that delays the rate at which WDs cool [lj. While 
observations of WD cooling in NGC 6791, a metal rich open cluster, shows that the 
WD luminosity function is well fit by evolutionary models only when including both 
sedimentation of neutron rich ions such as 22 Ne, which releases gravitational energy, 
and the release of latent heat of crystallization ji); [lj]. Modeling diffusion precisely 
in WDs allows one to estimate the rate of sedimentation and thus determine more 
accurately the age of a stellar cluster. 

While most stars will end as WDs, high-mass stars end as NSs. As a high mass 
star reaches the end of its lifespan its core can reach densities too high to be supported 
by electron degeneracy pressure. The electron Fermi energy increases to a point 
where inverse beta decay occurs and electrons are forced into the nuclei, 

p + e~ — > n + v. (1.1) 

It is well known that free neutrons are unstable and will beta decay with a half-life 
of just under 15 minutes. However, this reaction is not allowed to happen within a 
NS because the electron Fermi energy is too high. The neutrons in a NS are stable 
because of their environment jH . 

The typical size of a NS is approximately 12 km with a mass between 1.4-2.0 M Q . 
The maximum mass, size, and interior structure of a NS is very equation of state 
(EOS) dependent. The mechanism keeping a NS from further collapse is neutron 
degeneracy pressure and strong interactions. The outer crust of a neutron star is 
comparable in density to the interior of a WD and can have a complicated structure. 
If a NS continually accretes material from a companion star the infalling material 
can undergo nuclear reactions by rapidly capturing protons (the rp process). The 
material also increases the density of the liquid ocean layer until it crystallizes. The 
result is a crust of complex structure that is believed to have a top liquid ocean layer 
of low-Z nuclei with a bottom solid crust of high-Z nuclei. 

Phase separation of high- and low-Z ions in the liquid and solid crusts has been 
described and explicitly modeled with MD simulations by Horowitz et al. 2007 
j^j. Chemical separation between the liquid and solid crust is believed to change 
properties of the crust and impact many observables: thickness, shear modulus, and 
breaking strain of the crust. Possible observable effects include changes in the shape 
of a NS, radiation of gravitational waves, and properties of quasiperiodic oscillations 
observed in magnetar giant flares 0. 

Diffusion in Coulomb plasmas in the liquid phase has been well studied since 
the 1970s. Hansen et al. performed MD simulations for the one component plasma 



CHAPTER 1. INTRODUCTION 



3 



(OCP) which consists of ions interacting with pure Coulomb interactions and an 
inert neutralizing background charge density _|6|. More recent work has been done 
to understand diffusion in Coulomb crystals [7J. Particles exchanging lattice sites or 
the diffusion of imperfections in the crystal can be important mechanisms for the 
release of gravitational energy. In simulations, the diffusion in Coulomb crystals is 
dependent on the form of the inter-particle potential. The Lennard Jones potential 
with its hard-core contribution (r -12 ) tends to form a glass (sf with low diffusion 
rates while the classical Coulomb (1/r) potential allows ions to diffuse much more 
easily. 

Another commonly used interaction is that of the Yukawa potential. Diffusion 
for a Yukawa fluid has been simulated by Robbins et al. Q and Ohta et al. 10]. In 
a Yukawa fluid ions interact via a screened Coulomb potential Vij(r), 



2 



%(r) = ^-e~ r /\ (1.2) 

for two ions with charges Zj and Zj that are separated by a distance r. The OCP is 
equivalent to a Yukawa fluid, where all of the ions have the same charge (Zi = Zj), 
and the Thomas Fermi screening length A is very large. 

In WDs and NSs the motions of ions are typically considered to be classical be- 
cause of their large masses. However, at extreme densities, quantum corrections have 



been considered 111; Il2|; ll3j . In 2005, Daligault et al. studied quantum corrections in 
diffusion by modifying the interaction %(r) in a semiclassical approximation which 
includes the effects of zero point motion. It is the intent of this research to apply 
this semiclassical approximation to MD simulations of diffusion in crystals and to 
compare with work in which a purely classical approach has been used. 

In the interiors of WDs and the crusts of NSs the inter-ionic spacing, which is 
close to the size of the ion-sphere radius a = (3/Airn) 1 ^ 3 where n is the ion density, 
becomes comparable to the ionic thermal deBroglie wavelength A t h = \j2nh 2 / 'MT ', 



where M is the ionic mass [12| . As these values become comparable quantum effects 
become more important. 



We choose to define the parameter A = A th /\ / 2Tc 2 following ref. 11]. The 
dimensionless quantity A/a can be used as a measure of the importance of quantum 
effects in a system. For systems in which A/a C 1 interactions can be considered 
to be classical. For those where A/a ~ 1 quantum effects are extremely important. 
For WDs and NSs they lie between the two previous cases with < A/a < 1 and a 
semiclassical approximation to particle interactions can be used. The semiclassical 
addition to the ionic potential is given as a variation of the Yukawa potential and 
takes the form 

ZiZje 2 



r 



-r/X^_ e -r/Ay ( L3 ) 



This form for the interaction accounts for the extension of the particle by effectively 
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smearing out the ion on the scale of A. When A is small this form reduces to 
that of Eq. ( II. 2p . Such a form for the interactions makes the ions softer-core and 
should allow diffusion to more readily occur for larger values of A/a. Larger values 
of A/a can be achieved for lighter elements (eg. 4 He), colder temperatures, or higher 
densities. These effects are explored in more detail in chapter 4. 

In this paper we present our MD formalism in chapter [2l Diffusion constants 
and melting temperatures for classical and semiclassical simulations are presented in 
chapter |3j In chapter H] we discuss implications of this research on the understanding 
of structure within WDs and NSs. We conclude in chapter [5j 
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Chapter 2 
Formalism 



An MD simulation is a technique in which a system of interacting ions is evolved 
with time by integrating the ions' equations of motion. Given an initial set of posi- 
tions and velocities for all ions within the system, the time evolution is completely 
determined, in principle. Time evolution is based on Newton's law, F = ma, and the 
forces are obtained as the gradient of a chosen potential, which is a function of all the 
particle coordinates. Positions and velocities of ions are evolved through the velocity 
Verlet algorithm 15 L A detailed description of many aspects of MD simulations is 



given by Ercolessi [16 



2.1 Molecular Dynamics Simulations 

For the purposes of this research, we consider only a one component plasma 
(OCP). A star with near solar metallicity which has converted most of its original 
carbon, nitrogen, and oxygen into 22 Ne, might have of order 2% 22 Ne. The ratio 
of carbon to oxygen in the core depends on the rates for the 4 He(2a,7) 12 C and 
12 C(a,7) 16 reactions and is expected to be near one to one 17J. These conditions 
can lead to a WD with a core composition close to 49% 12 C, 49% 16 0, and 2% 22 Ne. 
Because of this high concentration of 16 0, we consider a OCP composed entirely of 
16 0. Later simulations may consider usin g a multi-component plasma to compare 
with previous work done by Hughto et al. [14J. 

Ions interact according to Eq. (jl.3p . The Thomas Fermi screening length A, 
for relativistic electrons, is A -1 = 2a 1 / 2 /^/V 1 / 2 where kp is the electron Fermi mo- 
mentum and is given as k F = (3ir 2 n e )V 3 and a is the fine structure constant. The 
electron density, n e , is equal to the ion charge density, n e = Zn, where n is the 
ion density and Z is the ion charge. For simplicity, we use the extreme relativisitic 
limit by neglecting the electron mass. However, this could become more important 
at lower densities, as this decreases A. 
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Simulations can be described by a Coulomb coupling parameter 

_ Z 2 e 2 



aT 



(2.1) 



Here T is the temperature of the system. T is the ratio of the Coulomb potential 
energy to the thermal kinetic energy, and it is known that a OCP crystallizes at 
r ?y 175. It is also known that the crystallization value of T may also depend slightly 



on the value of A 18 



One of the most fundamental timescales in plasma physics is that of the plasma 
frequency u p . Long wavelength fluctuations in the charge density can undergo oscil- 
lations at the plasma frequency. Therefore, in our simulations, time can be measured 
The plasma frequency depends on the ion charge Z and mass M, 



in units of cu p 1 



AnZ 2 e 2 n 
M 



1/2 



(2.2) 



All simulations are evolved with the velocity Verlet algorithm [15[ using time 
steps of At l/9u p . Periodic boundary conditions are used in each orthogonal 
direction. We do not use a cutoff distance for interaction and, therefore, calculate 
the interactions between all ions. The force on an ion is evaluated by summing over 
all other ions. In an effort to decrease finite size effects, all simulations use N = 8192 
with a box size L that is much larger than the electron screening length, L/2 = 8.92A. 
Temperatures are held nearly constant by periodically rescaling the velocities every 
2.36/u p (twenty time steps). Any bulk flow of ions is treated by calculating the 
center of mass velocity in all directions and subtracting from the velocity of each ion 
every 11% /uj p (one thousand time steps). 



2.2 Liquid Phase Simulations 

In liquids the diffusion constant D can be calculated from the velocity autocor- 
relation function, 

= Mto + t).^)) 
KJ (vM-Vjito)) 1 ; 

This averages over all ions j and over initial times to to minimize statistical errors. 
The velocity of the jth ion at time t is v, (t). The diffusion constant is calculated 
from the time integral of Z(t), 

D = - dtZ(t). (2.4) 
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With time Z(t) fluctuates about zero, and thus contributions to D are significantly 
reduced by a time t max = 240 /ui p . 

To ensure that simulations with the semiclassical approximation begin with sim- 
ilar initial conditions, a system with a classical interaction is allowed to equilibrate 
for a time of 1200/wp (10,000 time steps). The ions begin with random initial po- 
sitions at T = 150. Ions are given random initial velocities between ±3kT/M. As 
the simulation evolves the ions lose the random velocity distribution and assume 
one of a Gaussian distribution. The final positions and velocities for this system are 
written out to a file and used as the initial conditions for subsequent liquid phase 
simulations. All liquid phase simulations which include quantum corrections use the 
classical positions and velocities for their initial conditions. A value for A is then 
introduced and the systems are allowed to equilibrate for a time of 1200/u;p. Simu- 
lations are then to evolved for an additional 1200/w p during which the positions and 
velocities of ions are written to a trajectory file after every time step. Results for 
Z{t) and D(t) are given in section 1X21 



2.3 Solid Phase Simulations 

Diffusion in crystals is much smaller than in liquids making it difficult to calculate 
diffusion coefficients using Eq. (12. 3p . As Z{t) fluctuates about zero the integral in 
Eq. (12.41) involves sensitive cancellations of D{t) where Z(t) is positive and negative. 
A more effective method of calculating D{t) in solids is 

Bit) = (Mt + tJ-rWrc (2 . 5) 

where Tj(t) is the position of the jth ion at time t and the average is over all ions j 
and initial times to- The diffusion constant D is the large time limit of Eq. (12.51) . 

D = lim D(t). (2.6) 

t— >oo 

However, simulations are finite in time, and thus thermal oscillations about lattice 
sites can contribute to D(t) for small time limits of Eq. (12. 6p . An ion oscillating 
about a lattice site does not contribute to the net diffusion of the system, but it 
will have \rj(t) — (0) | nonzero and contribute to Eq. (12.61) . Therefore, thermal 
oscillations cause D(t) to differ from D for small t. It is convenient to define the 
quantity 

jy, t) = (e[|r,(^) -r,(t )| -^ c ]|r,(^) -r.Cto)! 2 ) ^ 

6t 

where t' = t + t . Eq. (12.71) introduces a cutoff radius R c for which ions are required 
to travel before contributing to the net diffusion. The cutoff radius is of order of the 
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lattice spacing and reduces contributions of thermal oscillations to D'(t). It has been 
observed by Hughto et al. 0] that D'{t) is approximately independent of t, even for 
moderate t, so that 

D « D\t). (2.8) 

It is also important to note that at arbitrarily large t, ions which diffuse distances 
nearly equivalent to the width of the box introduce error in Eq. (12.51) and (12. 6p due 
to periodic boundary conditions. However, diffusion is relatively slow in solids so 
this is often not a problem until very large t. 

Initial conditions are very important for determining D, as systems can possibly 
contain defects which may take a long time to equilibrate. Hughto et al. describe a 
simulation in [?J that suggests both WD and NS plasmas freeze into nearly perfect 
body-centerd cubic (bcc) crystals. Because of their result, we create initial conditions 
from a classically interacting system (A/a = 0) of a pure bcc crystal at a temperature 
equivalent to T = 300 and slowly warm it to V — 175. Finally, the quantum correc- 
tion is introduced to the interactions, and the system is allowed to equilibrate for 
1200/wp and then evolved for 4800/o;p (40,000 time steps) while the ions' positions 
and velocities are recorded every 23.Q/u p (200 time steps). 

Diffusion in a crystal at a temperature of T = 200 are also considered. The initial 
conditions begin from the same initial conditions as T = 175 simulations. However, 
when the T = 200 simulations are allowed to equilibrate the system is cooled off by 
adjusting the ions' velocities along with introducing the value for A. Systems are 
equilibrated for 1200/cUp and then evolved for 4800 /u p while the ions' positions and 
velocities are recorded every 23.6/u p . 

In both temperature regimes (r = 175, 200) Eq. (12. 7p is used to calculate the 
diffusion constant with t' ~ 4250u p . This allows the diffusion constant to be averaged 
over twenty configurations. The initial positions Tj(t ) of each configuration are 
separated by 23.6/u p . Final positions Tj(t') are separated from initial conditions by 
t m 4250/ajp. This is done to decrease statistical uncertainties in D. Histograms 
of number of ions diffused versus displacement are constructed and presented with 
diffusion results in chapter 13.31 



2.4 Melting Temperature 

Another aspect in which this research can provide insight is that of the melting 
temperature of a system that includes quantum corrections, equal to that of Eq. 
(11.31) . It is known that the melting temperature of a classically interacting system 



is T « 175 [18] . but little work has been done for systems which include quantum 
corrections. 

The equilibrium temperature is found by evolving a system of half liquid and half 
solid then visually checking if the system begins to favor one configuration. If the 
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system is out of equilibrium it will begin to liquify or crystallize accordingly, and 
the value of T is appropriately adjusted to raise or lower the temperature. Visual 
checks of the system are made using the program Visual Molecular Dynamics (VMD) 
created by the Theoretical and Computational Biophysics Group at the University 



of Illinois at Urbana- Champaign [23 



Initial conditions for both liquid and solid phases are created separately. Each 
phase consists of 8192 ions for a total of N = 16384 in the combined system. The 
liquid phase is created the same way as that described in section 12.21 The system 
is then cooled to T = 175 and equilibrated for 1200/cjp. The solid phase is acquired 
from the same initial conditions described in section 12.31 at F = 175 as this con- 
figuration was readily available from previous simulations. The separate halves are 
brought together by adding a value of L/2 — 8.92A to the Z-component of the liquid 
configuration. Periodic boundary conditions along the Z-axis are altered to allow for 
a now rectangular box. A sample configuration of a two-phase system is given in 
Figure I2TT1 

When the two halves are brought together the boundary between phases will 
have ions which are closer than would typically occur. This is because the two 
halves are equilibrated separately. Combined systems require time to equilibrate as 
the particles near boundaries adjust to their new neighbors. Shown in Figure 12.21 
is the range of temperatures and equilibration time given to bring the classically 
interacting two-phase simulation near the melting temperature. 
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Figure 2.1: Sample configuration of N=16384 ions showing a liquid phase on the left 
and solid phase on the right. Figure prepared with VMD 
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Figure 2.2: Range of temperatures used to bring the two-phase mixture of classically 
interacting (A/a = 0) 16 ions near to its melting temperature. The system was 
given a total simulation time of t ~ 6100Av p (52,000 time steps). 



In order to understand how quantum corrections affect melting temperature, two- 
phase system with A/a = 0.335 and A/a = 0.502 were evolved. Initial conditions for 
simulations which include quantum corrections are created using the same method 
for initial conditions as the classical simulations described in the previous paragraph. 
Both liquid and solid phase configurations are equilibrated separately for 1200/w p as 
before, except now, during the equilibration phase Eq. fll.3p is used for the inter-ion 
potential. The separate halves are also equilibrated at cooler temperatures than the 
classical simulation. A larger T is chosen for this step because quantum corrections 
can melt the system due to zero point motion. This will be discussed in greater 
detail in chapters [3] and HI The two phases are then brought together as before and 
allowed to equilibrate while manually adjusting T to keep the system half solid and 
half liquid. 

Both classical and semiclassical simulations are brought near to their melting 
temperatures by visually checking if the systems are changing phase and adjusting 
the value of T. Once the two-phase system is near the melting temperature it is 
then evolved microcanonically, at a constant energy, by no longer rescaling the ions' 
velocities. This method allows the release or absorption of latent heat allowing 
the system to adjust its own temperature and liquify or crystallize as needed. The 
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interface between the liquid and solid, along with the temperature, will exponentially 



reach an equilibrium state [16|. Once the system has self-equilibrated, the melting 
temperature is inferred from the kinetic energy of the system. Results for these 
simulations are given in section 14.21 
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Chapter 3 



Results 



In this chapter we present results for all Molecular Dynamics (MD) simulations. 
Typical computation time for N = 8192 ions and a simulation time of t ~ 4800/a;p 
(40,000 steps) is approximately two-and-a-half days. Trajectory file sizes for liquid 
simulations are typically 8.3 GB in size. The authors wish to thank Indiana Univer- 
sity's Scholarly Data Archive (SDA), formerly known as MDSS, for providing storage 
of data files. Calculation of the velocity autocorrelation function, which is compli- 
cated and very memory intensive, required on average thirteen hours. All simulations 

TM 

were performed on a Dell Optiplex 760 desktop. This 64-bit machine contained 

TM 

an Intel® Core 2 Duo CPU E8400 with 3.0 GHz processor speed along with 4.0 
GB of RAM. 

3.1 Potential Energy per Particle 

An important test for MD simulations is to evaluate the energy involved within 
a system. The average potential energy per particle (V) is calculated using 



where Vij is equal to Eq. (JT72J) for classical simulations or Eq. ( II. 3p for the semiclassi- 
cal simulations. Figure I3TT1 and Table I3TT1 show the change in (V) as a function of A/a 
for the three temperatures used throughout this research. In crystals near the melt- 
ing temperature thermal oscillations are large and ions come closer together and lead 
to higher total potentials. Thermal oscillations are smaller for lower temperatures- 
larger T-and ions stay farther apart. This decreases the total potential energy of the 
system. While the colder systems begin with a lower potential energy per particle, 
the characteristic to note is that the potential per particle decreases for increasing 



A 1 




(3.1) 



i<3 
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Table 3.1. (V) versus A/a for All Systems 



A (fm) 


A/a 


T = 150 


T = 175 


r = 200 






(MeV) 


(MeV) 


(MeV) 








26.816 


26.781 


26.773 


1 


0.067 


26.816 


26.781 


26.773 


2.5 


0.167 


26.816 


26.781 


26.773 


4 


0.268 


26.806 


26.773 


26.764 


5 


0.335 


26.779 


26.749 


26.738 


6 


0.402 


26.723 


26.713 


26.685 


7.5 


0.502 


26.581 


26.572 


26.564 


10 


0.670 


26.202 


26.194 


26.187 



A/a. The average decrease in energy for all simulation temperatures is 0.596 MeV. 
Both solid systems change phase due to zero point motion for large A/a, but the 
system's energy per particle changes very little. As the quantum influences become 
larger the potential becomes more attractive at short distances which decreases the 
potential per particle. 
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Figure 3.1: Potential energy per particle, (V), versus A/a for all simulations. As 
quantum effects become larger the energy per particle continues to decrease due to 
decreased inter-ion potential. 
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3.2 Diffusion in Liquids 

Figure I3.2al shows the velocity autocorrelation function Z(t), Eq. ( 12.31) . for a 
classically interacting system and several intermediate values of A/a. Smaller values 
of A/a are nearly indistinguishable from that of A/a = and are not shown. This 
figure shows that for increasing A/a, Z(t) tends to 'lag' behind that of a classically 
interacting system. The velocity autocorrelation function oscillates with a frequency 
near u p , but zero point motion decreases the frequency with which Z(t) oscillates. 
Figure I3.2bl shows how the frequency of the velocity autocorrelation function de- 
creases as A/a is increased. The plasma frequency, Eq. (12. 2p . is used to scale all 
frequencies. 
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Figure 3.3: Integral of Z(t). (l2.4p . for several T = 150 liquid simulations. Only values 
of A/a which dramatically increased the diffusion are represented. 



Figure 1331 shows the integral of Z(t), Eq. (12.41) . Values for D are given in units 
of id p a 2 . Large fluctuations in D(t) can be seen for small time t due to large changes 
in Z(t). As t increases to t max D(t) converges to that of the true value. 

Diffusion constants for a one component plasma (OCP) in the liquid phase can 
be scaled by 



which is given by Hansen et aVs fit to their original MD results for diffusion [6[. 
Figure l3~4l and Table [3721 show the scaled diffusion constants as a function of A/a. 
Diffusion in simulations with values of A/a < 0.3 deviates very little from Hansen 
et a/.'s classical, theoretical diffusion given by Eq. (13. 2(1 . As quantum corrections 
increase, the diffusion slowly increases to twice the value predicted by Hansen et al. 
when A is two-thirds of the ion-sphere radius. The value of A/a does not appear to 
be extremely important to the diffusion in a liquid system. The ions are allowed to 
move freely around one another in a liquid so the importance of zero point motion is 
reduced. Large values of A/a help increase the diffusion, but the diffusion in a liquid 
is already relatively large. 
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Table 3.2. Diffusion Constants for Liquid Phase Simulations 



A (fm) 


A/a 


D/( 


u t 


a 2 ) 




D/D 








3.632 


X 


10" 


-3 


0.965 


1 


0.067 


3.564 


X 


io- 


-3 


0.947 


2.5 


0.167 


3.539 


X 


io- 


-3 


0.940 


4 


0.268 


3.742 


X 


io- 


-3 


0.994 


5 


0.335 


4.144 


X 


io- 


-3 


1.101 


6 


0.402 


4.654 


X 


io- 


-3 


1.237 


7.5 


0.502 


5.677 


X 


io- 


-3 


1.508 


10 


0.670 


7.417 


X 


io- 


-3 


1.970 



2.5 



2 



1.5 

o 

Q 

\ 

1 



0.5 





0.2 0.4 0.6 0.8 1 

A/a 

Figure 3.4: Diffusion for liquid simulations scaled by D , Eq. (13. 2p . 
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3.3 Diffusion in Solids 

We begin by showing several histograms of displacements \rj(t+to)—Tj(to)\ at the 
end of the simulations, t m 4800/w p . These are computed by counting the number of 
ions that have moved a given distance in a time t. Figure 13.51 is the histogram for a 
purely classical system and shows a large central peak at small distances, which is due 
to each ion remaining at its original lattice site. The width of this peak corresponds 
to thermal oscillations of ions about their respective lattice sites. The smaller peak 
centered at 1.8a corresponds to ions that have "hopped" to neighboring lattice sites. 
A description of how ions move within a crystal lattice is given by Hughto et at 
(2011) [7J. In this simulation a total of 372 ±19 ions (~ 4.5% of the total number of 
ions) have moved farther than 1.07a from their original lattice site. This distance is 
also the cutoff distance used in Eq. (12. 7p . 



1000 



2; 




2 4 6 

|r(t)-r(0)|/a 

Figure 3.5: Histogram of displacements \rj(t + to) — rj(to) \ in units of the ion-sphere 
radius a for a classically interacting crystal lattice after a time t ~ 4800/u;p. 

Figure l3.6al shows the displacement histogram for A/a = 0.335. The central peak 
for this histogram is smaller than the classical case, and the second peak is much 
larger than before. Many more ions are allowed to hop not just to neighboring lattice 
sites but to much farther distances. The distended peak out to ~ 6a shows that 
many more ions have diffused in this semiclassical case. In this simulation 4405 ±66 
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ions (~ 54%) have diffused farther than the cutoff distance. As A/a increases, the 
inter-ion potential decreases at intermediate distance and becomes more attractive 
at short distances allowing ions to diffuse more easily. As A/a increases even further, 
the system can then be melted by zero point motion. This is shown in Figure l3.6bl as 
there is only one extended peak centered around a distance of ~ 8a. At this point the 
zero point motion of ions has melted the crystal structure and ions can flow freely. 
Zero point motion is needed for ions to diffuse around neighboring particles, and 
large quantum corrections decrease the potential between ions, allowing diffusion to 
more readily occur. 





2 4 6 10 20 

|r(t)-r(0)|/a |r(t)-r(0)|/a 
(a) Histogram of displacements for A/a = 0.335. (b) Histogram of displacements for A/a = 0.402. 



Figure 3.6: Histogram of Displacements for the Semiclassical Potential 



Figure 13.71 and Table 13.31 show the diffusion constants for all solid simulations. 
Diffusion in solids cannot be scaled by Eq. (13. 2p because this equation only predicts 
diffusion in the liquid phase. For simulations with T = 175 the diffusion changes 
very little while A/a < 0.3. At A/a m 0.3 the diffusion increases by a factor of 20 
from a purely classical system, which can be seen as a large upturn in Figure 13.71 
For A/a > 0.4 the crystal structure is completely melted and the diffusion increases 
by a factor of 500. Once the system has melted, increasing A/a further does not 
have a significant effect and the diffusion only increases slightly. 

Crystal simulations with T = 200 exhibit similar characteristics. However, the 
system requires a larger value of A/a before melting occurs. For a system to melt 
a combination of thermal and zero point motion is required. For the colder system 
the thermal motion is reduced; therefore, more quantum motion is required before 
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melting can occur. This explains why A/a must be larger to melt the colder system. 
We recognize the idea that the system may be super heating and address this later 
by calculating the melting temperature directly. 

Of particular note, the simulation with T = 200 and A/a = 0.167 exhibited very 
little diffusion. During the simulation several ions drifted farther than the cutoff 
distance and were considered to have diffused. Each of these ions subsequently 
returned to their original, vacant lattice sites within several more time steps. At the 
time of the simulation's end two ions had moved farther than the cutoff distance, 
giving a value for the diffusion that is two orders of magnitude smaller than any 
other simulation, D/u p a 2 = 1.6 x 10~ 8 . It is likely that these two ions would have 
returned to their original lattice sites given more time. 
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Figure 3.7: Diffusion constant D versus A/a for solid phase simulations. 




CHAPTER 3. RESULTS 



20 



Table 3.3. Diffusion Constants for Solid Phase Simulations 
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3.4 Melting Temperature 

Multi-phase systems are brought close to their equilibrium temperature by man- 
ually changing the value for T. Systems are then evolved at constant energy allowing 
the system to find an equilibrium temperature without external influences. When 
running simulations at constant energy a calculated value for T is written to a file 
every 10 steps. Averages of T are then made for every 1000 simulation steps and plot- 
ted versus time t. This shows how the system evolves and changes the temperature 
to reach an equilibrium point. Figure 13.81 shows temperature changes throughout 
the simulation for a classical system while running at constant energy. Simulations 
are considered to have reached the melting temperature when changes in T become 
small. 
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Figure 3.8: Average T versus simulation time in a classical simulation. Each data 
point is the average of 100 values of T spaced ten time steps apart. 

In our classically interacting system (A/a = 0) the temperature initially behaves 
as expected, and increases towards the melting temperature. However, after a time 
t = 2100/cUp (18,000 time steps) the temperature begins to fluctuate. This system 
was given considerable time for fluctuations to diminish. To contest this the melting 
temperature is taken as the long-time average of many values. Only the final t ~ 
4000/ojp (34,000 time steps) are used in the average of T. Our inferred melting 
temperature of the classically interacting system is T = 187.6 ±1. The uncertainty 
is a conservative estimate that includes both the maximum and minimum of the 



fluctuations. This temperature is colder than the expected value of T 175 [18 
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Differences in our melting temperature from that of the expected value could be due 
to finite size effects and screening effects. 
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Figure 3.9: Average T versus simulation time in a semiclassical simulation with 
A/a = 0.335. Each data point is the average of 100 values of T spaced ten time steps 
apart. 

To understand how quantum corrections affect the melting temperature, a multi- 
phase simulation with A/a = 0.335 has also been performed. This value was chosen 
because it is the stage in the solid simulations where diffusion was beginning to oc- 
cur on a large scale basis. Figure 13.91 shows the temperature of our A/a = 0.335 
simulation versus time. Fortunately, this simulation does not have the appreciable 
temperature fluctuations that the classical system shows in Figure 13. 8[ and, there- 
fore, was run for a fraction of the computing time. From the A/a = 0.335 simulation 
we infer a melting temperature of T m = 206.5 ±1. A 10% increase in T is needed 
to bring the system to an equilibrium condition. Again, the uncertainty is a con- 
servative estimate which considers the maximum and minimum of the temperature 
fluctuations. 
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Table 3.4. Melting Temperature versus Quantum Corrections 



A (fm) 


A/a 


r 

L m 





0.0 


187.6 ±1 


5 


0.335 


206.5 ±1 


7.5 


0.502 


253.2 ±1 



254 



252 



250 



248 



246 



500 1000 1500 2000 

t*CJ 

p 

Figure 3.10: Average T versus simulation time in a semiclassical simulation with 
A/a = 0.502. Each data point is the average of 100 values of T spaced ten time steps 
apart. 

Finally, a multi-phase system with A/a = 0.502 has also been run. This value was 
chosen because, at this stage, the zero point motion of the ions is large enough to melt 
both of our solid simulations. Shown in Figure 13". 101 is the temperature change of the 
A/a = 0.502 simulation while evolving at constant energy. The melting temperature 
we infer from this data is r m = 253.2 ±1. This corresponds to an increase of 
35% in the melting temperature. Table [3741 gives the inferred melting temperatures 
versus quantum corrections. The implications of quantum corrections to the melting 
temperature are discussed further in chapter 14.21 
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Chapter 4 
Discussion 



4.1 Quantum Dependence on M, T, and n 



In this section we go beyond the results of our Molecular Dynamics (MD) simu- 
lations and discuss the implications that this research has on our understanding of 
white dwarfs (WDs) and neutron stars (NSs). Recall that A=A^/ V2tt 2 where A#, 
is the ionic thermal deBroglie wavelength and is given by A th = ^2ixh 2 /MT. Also 
recall that the ion-sphere radius is found by a = (3/47m) 1 / 3 . By combining these 
with a form of Eq. (12. ip that has been solved for T, we find 



A 



h 2 



r 



71 MZiZjC 2 



An 



1/3- 



1/2 



n 



1/6 



(4.1) 



where M is the ionic mass, Zi the ionic charge, T the Coulomb coupling parameter, 
and n the ion density. From our simulations, specific values of A/a have a respective 
associated diffusion, and in this chapter we consider several dependences of Eq. (14.1 ft 
to understand how each one can affect the diffusion of a system. 



4.1.1 Composition 

Here we discuss the A/a dependence on the composition of a system. Various 
compositions will change M, Zi and Zj in Eq. ( 14. ip . This equation shows that A/a 
oc (ZjZjM) -1 / 2 . Therefore, heavier ions have smaller thermal deBroglie wavelengths 
compared to the ion-sphere radius. Calculations are done with T = 175 and a 
relatively high ion density of n = 7.18 x 10 -5 fm~ 3 . This density corresponds to a 
mass density of 1.91 x 10 12 g cm -3 and was chosen for historical reasons. Ionic mass 
M is chosen using common stable isotopes of each element. Figure 14.11 shows how 
A/a decreases with increasing atomic mass and ionic charge. 

It can be seen in Figure 14.11 that light elements can have large ionic thermal 
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deBroglie wavelengths compared to the ion-sphere radius. The ion 4 He can have a 
thermal deBroglie wavelength more than twice the ion-sphere radius and 1 H can have 
a wavelength ten times the ion-sphere radius. These light elements are in the extreme 
quantum regime where small quantum corrections to the potential, such as that used 
in Eq. (11.31) . may be insufficient. On the other hand, heavy elements have smaller 
thermal deBroglie wavelengths. Ions such as 56 Fe, at this density, have an associated 
A/a ~ 0.05 which puts it in the classical regime. It is the intermediate mass ions 
(eg. 12 C and 16 0) which fall into the semiclassical regime, where quantum corrections 
begin to become important, and with which this research is most concerned. 
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Figure 4.1: Values of A/a versus atomic number (Z, M). The bottom graph con- 
tains the same data as the top graph, but it is scaled to the range < A/a < 1. 
Calculations are performed at T = 175 and an ion density n = 7.18 x 10~ 5 fm -3 . 



4.1.2 Temperature 

Next we consider the dependence of A/a on temperature, for which T is strongly 
correlated through Eq. (12. ip . Increasing T is analogous to decreasing the temperature 
and is very similar to a cooling WD. The density and composition of the system, 
aside from sedimentation, should change very little over time, but as it cools the value 
for T will continue to increase. Eq. (14. ip shows that A/a oc T 1 / 2 . In calculations, 
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we evaluate Eq. (14. ip using 16 O and an ion density of n = 7.18 x 10~ 5 fm~ 3 . Figure 
14.21 shows that at T = 150 we find A/a k 0.3, which is where quantum effects 
began increasing diffusion in our liquid simulations. This implies that a liquid may 
experience increasing diffusion as it cools and nears crystallization. This could also 
further chemical separation in the liquid crust of a NS 0; [l^. Furthermore, at 
T = 175, A/a is approximately 0.32, which is where our simulation's diffusion had 
increased by a factor of 20. Finally, our T = 175 simulations experience a complete 
phase transition at A/a > 0.4. Eq. (14. ip does not evaluate to 0.4 until T « 265. 
This semiclassical approximation implies that a WD could reliquify as it cools. 
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Figure 4.2: Values of A/a versus temperature T. Calculations performed with 16 
at an ion density n = 7.18 x 10~ 5 fm~ 3 . 



4.1.3 Density 

We now consider the A/a dependence on ion density n. This is used to understand 
how A/a behaves as you move to the interior of a WD or deeper into the crust of a NS. 
We first discuss the range of densities found in a WD or NS. WD stars typically have 
densities of p ~ 10 6 g cm -3 , but in an attempt to characterize ignition conditions 
of type la supernovae Lesaffre et al. 2(| find that for a range of stellar masses the 
central density for ignition for WDs falls within 2 x 10 9 — 5 x 10 9 g cm -3 . A WD 
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should not be expected to reach central densities above this or it would otherwise 



become a supernova [20J. However, a NS is comprised of densities far beyond that of 
a WD. In a system where a NS is accreting mass from a companion star it will have a 
top liquid ocean layer; an outer crust which consists of a 1 H burning layer at a density 
p ~ 10 5 g cm -3 ; a level of 4 He at p ~ 10 6 — 10 8 g cm -3 which can burn unstably 
because it is strongly degenerate; and a bottom layer with p ~ 10 11 g cm -3 where 



neutronization becomes relevant, ions are neutron rich, and neutron drip begins [21 
The inner crust densities reach p ~ 10 13 — 10 14 g cm -3 and become a homogeneous 



mixture of n, p, and e with few percent protons at the transition to the core [22 



The bottom of the inner crust is near normal nuclear density p = 2.7 x 10 14 g cm 3 
21 . 

It is also important to mention the densities at which pycnonuclear reactions oc- 
cur {ppyc) and densities for neutronization (p n ). Pycnonuclear reactions are heavily 
dependent on density rather than temperature. Thermal vibrations of ions about 
their lattice, along with the probability to tunnel through the repulsive Coulomb 
barrier, can lead to nuclear reactions. Reactions set in quickly at p > p pyc . Typical 
densities for pycnonuclear reactions are p pyc ~ 10 6 , 10 9 , and 10 10 g cm -3 for burning 
1 H, 4 He, and 12 C respectively. Because of their larger Coulomb barriers heavy ele- 
ments require higher p pyc . For each ion there is also a threshold p n of the density for 
which neutronization occurs. For 1 H, 4 He, and 12 C this density is p n ~ 1.2 x 10 7 , 
1.4 x 10 11 , and 3.9 x 10 10 g cm -3 respectively. For all ions mentioned, p pyc <C p n and 
pycnonuclear reactions occur before neutronization 

The thermal deBroglie wavelength, = ^J2nh 2 / 'MT ', does not explicitly depend 
on the density of the system. However, the temperature T does have a density 
dependence. Eq. 04. ip shows that A/a oc n 1 / 6 . The size of the thermal deBroglie 
wavelength is decreased for dense systems, but not as rapidly as the ion-sphere radius. 
For large densities the ratio A/a continues to increase and quantum corrections 
become more important. 
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Figure 4.3: Values of A/a versus density p. Calculations are performed with 16 0. 



To get an initial perspective on what density range quantum effects become im- 
portant to diffusion we use Eq. (14. II) and solve for the ion density n. For a given ion, 
16 0, at several T, the density dependence of A/a is plotted in Figure 14731 This figure 
shows a wide range of densities needed for A/a to increase from the classical regime, 
through the semiclassical, and into the quantum regime. An interesting thing to note 
in Figure 14751 is that for liquids, T = 150, higher densities are consistently needed to 
achieve the same A/a as crystal structures because of the higher temperatures. 

From our simulations, we can go further and now associate diffusion with density. 
Figures 14741 and 14751 are the same as Fi gur es [3741 and 13771 resp ect ively and show diffusion 
versus A/a for all simulation temperatures. However, Eq. ( 14. ip has been used to 
calculate the density at each value of A/a, at its respective T. These densities are 
also given in Table 13.31 All simulations were composed of a OCP of 16 O, therefore, 
all calculations for densities in these figures and table use a OCP of 16 as well. 
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Figure 4.4: Diffusion constants for liquid phase simulations with densities calculated 
from Eq. (14. ip . Calculations performed with 16 and T = 150. All densities are 
given in units of g cm -3 . 
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Table 4.1. Calculated Densities from Eq. ( 14. II) . 
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Figure 4.5: Diffusion constants for solid phase simulations with densities calculated 
from Eq. (14.11) . Calculations performed with 16 0, T = 175, and T = 200. All 
densities are given in units of g cm -3 . 



Finally we describe the diffusion dependence not only on density but on both 
density and composition. It has already been shown that light elements can have 
relatively large thermal deBroglie wavelengths, and that as the density of a system 
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increases quantum corrections become more important allowing diffusion to occur 
more readily. We do this to understand what densities various ions require before 
quantum effects become important. This work is complementary since sedimentation 
in WDs [l| and phase separation in the crust of a NS [^j are expected to occur for 
different compositions. 

This is done by choosing a specific ion (Z, M), holding T constant, holding A/a 
constant, and then solving for the density n in Eq. (14. ip . Figure 14.61 shows the 
diffusion dependence on density for common elements within a WD or NS crust (eg. 
4 He, 12 C, 16 0, and 56 Fe). 

It can be seen that light ions require low densities before quantum corrections 
become important. For example, 4 He at a V = 175 requires a density p ~ 2 x 10 6 
g cm -3 before zero point motion can melt the system. This density occurs in both 
WDs and NS crusts and is below the threshold densities p pyc and p n . 

Figure 14.61 shows that 56 Fe requires p ~ 10 17 g cm -3 for all simulation temper- 
atures before zero point motion becomes important. However, 56 Fe has a density 
for neutronization (p n m 1.14 x 10 9 g cm" 3 ) that is below its pycnonuclear reaction 
density jij. Therefore, neutronization of 56 Fe occurs before nuclear reactions occur. 
The density needed for quantum corrections to become important is far beyond the 
neutronization density and the density achieved by a NS. 
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Figure 4.6: Diffusion constants for all simulations as a function of density. The type 
of line (i.e. solid versus dashed) represents the temperature of the simulation. The 
color of the line represents the ion (Z, M) used in the density calculations. 



Intermediate mass ions 



i.e. 



12/ 



C and O) need much lower densities than Fe 



before quantum corrections become important. A solid system composed mostly of 
12 C near the melting temperature, T = 175, could remain a liquid if the density is 
over ~ 1 x 10 11 g cm" 3 . Such a density is beyond that of a WD, but it can be found 
in the bottom of the outer crust of a NS. However, this density is greater than p pyc 
for 12 C and nuclear reactions should occur first. We find that for both 12 C and 16 
the densities required for relevant quantum corrections are beyond that of the WD 
peak central densities stated by Lesaffre et. al 20|. However, zero point motion is 



relevant in the bottom of the outer crust of a NS where such densities can be found. 
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4.2 Quantum Effects on Melting Temperature 

We now discuss how quantum corrections influence the melting temperature. 
Figure 14.71 depicts the melting Coulomb parameter T m of our two-phase simulations 
versus A /a. 




A/a 

Figure 4.7: Melting Coulomb parameter T m for two-phase simulations. All simula- 
tions contain N = 16384 ions. Error bars are ±1 in magnitude and are a conservative 
estimate on the statistical error on T m . 



We find that T m increases rapidly for A/a > 0.3. This trend in Figure I4T71 may be 
much smoother than represented if simulations were to be run at more values of A/a. 
Modifying the melting temperature of material in a WD could play an important 
role in its evolution because the cooling mechanism depends on the phase of matter. 
This alters when the radiation of the latent heat of fusion will delay cooling 24 . 
However, quantum corrections do not impact melting temperature significantly until 
A/a ~ 0.335 which corresponds to a density of 2.3 x 10 12 g cm~ 3 for 16 at V = 175. 
This is higher than the densities reached in a WD, and, therefore, the ions in C and 
O WDs can be considered classical. For light elements, such as 4 He, the density 
required is only of order 10 6 g cm -3 , and melting temperatures could be altered by 
quantum effects. This implies that a He white dwarf would not freeze, and there 
would be no delay in the cooling curve as it radiates latent heat. 



CHAPTER 4. DISCUSSION 



34 



Table 4.2. Values of r s Corresponding to A/a 



A (fm) 


A/a 


T = 175 


T = 206.5 


T = 251.5 


1 


0.067 


12400 


14640 


17830 


2.5 


0.167 


1990 


2340 


2850 


4 


0.268 


780 


920 


1110 


5 


0.335 


500 


590 


710 


6 


0.402 


350 


410 


500 


7.5 


0.502 


220 


260 


320 


10 


0.670 


120 


150 


180 



We now compare our melting temperature results to that of Jones and Ceperley 
(1996) [25[. Their work used path integral Monte Carlo (PIMC) simulations to study 
the OCP at finite temperature and directly calculate the importance of quantum 
effects on melting temperatures in two-phase systems. Simulations are described by 
the dimensionless ratio r s = a/ao, where a is the ion-sphere radius and ao is a natural 
length scale given by ao = H 2 /mZ 2 e 2 . Large values of r s correspond to more classical 
systems while smaller values are more quantum in nature. In Jones and Ceperley's 
Fig. 1, they find that the melting temperature of a OCP is similar to the classical 
prediction for a large portion of their phase diagram. This figure also shows that 
for a given temperature as the density increases a OCP will solidify. Increasing the 
density even further increases quantum effects and can reliquify the system. This is 
in agreement with our results in chapter 14.1.31 At very high densities A/a becomes 
large and the melting temperature becomes very low, T m approaches infinity. Even 
at T = the system can remain a quantum fluid. 

Given in Table H~2l is a comparison of r s to A/a for all simulation temperatures. 
Note that large values of A/a correspond to small values of r s . However, from our 
simulations when A/a = 0.335 we find a 10% increase in T m from the classical case, 
see Table 13.41 At the corresponding density from chapter 14.1.31 this gives a value of 
r s ~ 590. Jones and Ceperley do not find a noticeable change in T m until somewhat 
smaller values, r s ~ 200. In the case of A/a = 0.502 we find a 35% increase in T m . 
For this condition's corresponding density we find r s ~ 320. Comparing our data 
to that of Fig. 1 in Jones and Ceperley our predictions for the melting temperature 



fall between the semiclassical prediction of ref. |26| and that of the full quantum 
calculations. 
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Chapter 5 

Summary and Conclusions 



Quantum corrections to the inter-ion potential can be important for dense as- 
trophysical objects and can affect the amount of diffusion in these systems. We 
have performed MD simulations of a OCP to understand the impact of quantum 
corrections. Diffusion coefficients for liquid and solid phase simulations have been 
calculated. Quantum corrections in the liquid configurations do not dramatically 
affect diffusion. Quantum corrections depend on the parameter A which is related to 
the ionic thermal deBroglie wavelength of an ion. Simulations show that the diffu- 
sion coefficient only increases by a factor of two for A/a = 0.670. However, quantum 
corrections in solid systems are much more important. We find that for systems near 
the melting temperature increasing A/a to 0.335 increases the diffusion coefficient 
by a factor of 20, and for A/a = 0.402 the system can be melted. 

Shortly after crystallization in the core of WDs, where T ~ 200, quantum effects 
are unimportant unless p > 1.5 x 10 12 g cm" 3 for 16 0, 1.5 x 10 6 for 4 He, and 2.7 x 10 17 
for 56 Fe. Quantum effects for ions between 16 and 56 Fe are likely to be small at 
WD densities. Quantum corrections for 4 He could be larger if 4 He survives to high 
densities. In NSs, quantum effects for 16 should be considered as the essential 
density can occur in the inner crust. However, the density needed for zero point 
motion to be significant for 56 Fe is well beyond that reached by WDs or NSs. 

It is important to understand the role of quantum corrections to the melting 
temperature in dense systems. Altering the melting temperature of a WD determines 
how long the star has to cool before it crystallizes, and in the crust of a NS it 
can affect the structure and depth at which crystallization occurs. We determined 
melting temperatures by performing two phase MD simulations where both liquid and 
solid phases are equilibrated simultaneously. We find that increasing the quantum 
corrections to increases r m . For A/a = 0.335 we find Y m = 206.5 ±1, a 10% 
increase, and for A/a = 0.502 we find T m = 250 ±1, a 35% increase. For 16 this 
corresponds to densities of 1.4 x 10 12 and 8.7 x 10 12 g cm -3 respectively, which can 
be found in NSs. For larger A/a the value of T m may become very large. Indeed, 
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at very high density there can be a quantum fluid that remains liquid even at zero 
temperature. 

In conclusion, we find quantum corrections to be small in WDs unless considering 
light elements such as 4 He. The quantum effects in the crust of NSs should be 
considered for light-to-intermediate elements (e.g. 12 C and 16 0), if these ions survive 
to high densities in the inner crust. For ions heavier than 56 Fe the densities required 
for quantum corrections appear to be higher than achieved by a NS and the ions can 
be considered classical. 
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